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ABSTRACT 



In this thesis, the use of the Wentzel-Kramers-Brillouin (WKB) Theory to obtain the 
solution to the Helmholtz Equation governing the acoustic normal modes is examined. 
Specifically, uniformly valid WKB solutions for four classes of acoustic normal modes 
in the ocean are derived and the accuracy of the WKB approximation is tested against 
some exact solutions. It is found that this inherently high frequency technique has an 
appreciable accuracy even at a frequency of 1 Hz. A product of this thesis is a computer 
program that solves for the WKB modes for an arbitrary sound speed profile. 




THESIS DISCLAIMER 

The reader is cautioned that computer programs developed in this research may not 
have been exercised for all cases of interest. While every effort has been made, within 
the time available, to ensure that the programs are free of computational and logic er- 
rors, they cannot be considered validated. Any application of these programs without 
additional verification is at the risk of the user. 
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I. INTRODUCTION 



A. BACKGROUND 

To model sound propagation in the ocean we can use several theories, namely Ray 
Theory, Parabolic Equation Approximation and Normal Modes. The Ray Theory sol- 
ution is an asymptotic geometric-optics solution, obtainable using simple ray tracing 
techniques. It provides a simple physical description on how sound is transmitted 
underwater. However, it neglects sound diffraction and thus needs corrections near 
caustics and turning points. Such corrections can be very complicated mathematically. 
The Parabolic Equation method is less physical but is a full-wave solution. An asset of 
Parabolic Equation is its capability to handle variable bottom bathymetry well. Its lim- 
itation is that it does not accurately model sound energy propagating in steep angles. 
Normal Mode Theory gives a physical full-wave solution to the wave equation. It has 
some computational difficulties in handling range varying sound speed fields. 

A computational difficulty associated with normal mode models applied to a range- 
dependent ocean is that the normal modes associated with a great number of profiles 
must be computed and the results stored. Large storage and processing time are re- 
quired if straight-forward numerical methods such as finite differences are used to com- 
pute the normal modes (Chiu and Ehret, 1990). The use of finite differences methods 
requires the discretized mode functions and their eigenvalues at a great number of points 
to be stored. Moreover, at high frequencies these methods require the computation of 
eigenvectors and eigenvalues of large matrices, which can result in significant increases 
of processing time and numerical noise in the solution. In this thesis we examine an 
asymptotic expansion method that has the potential to overcome this difficulty. The 
method is called Wentzel-Kramers-Brillouin (WKB) theory. 

B. THE NORMAL MODE APPROACH 

The wave equation governing the sound pressure p in the ocean is 



where c = c(z; r, 6) is the sound speed and z is the vertical coordinate, r the range and 6 
the azimuthal angle. In the three-dimensional coupled mode model of Chiu and Ehret 
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(1990). the pressure is expressed as a linear combination of the local normal modes Z„ 
such that 



oo 



p{r,e,z) = Y,^,{z-,r,e)P,{r,e)T{i) . 



For a harmonic frequency time dependence 




where &> is the acoustic angular frequency, the local modes at each horizontal location 
(r, 6) are required to satisfy 



and the appropriate boundary conditions. This equation is usually known as the 
Helmlioltz Equation. The constant k„ is the horizontal component of the wavenumber 
vector whose magnitude is given by 




In general, there are many possible values of k„ (eigenvalues) satisfying the Helmlioltz 
Equation. For each k„ there is an associated mode Z„ (eigenfunction). In a range- 
dependent sound speed field this eigenvalue-eigenfunction problem must be solved at a 
great number of horizontal grid points. 

C. THESIS OBJECTIVES AND OUTLINE 

The main objective of this thesis is to examine the use of the WKB theory to solve 
the Helmholtz Equation governing the acoustic normal modes. This includes: 

1. the development of the various WKB formulae for the four classes of normal modes 
which can exist in single duct/channel environments, 

2. the parameterization of each class of normal modes using a minimum number of 
parameters, 

3. the quantification of errors in the WKB solution through comparisons to three 
exact solutions. 

A product coming out of this thesis is a computer program solving for the mode 
parameters for an arbitrary sound speed profile. The incorporation of this code in the 
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three-dimensional coupled mode model of Chiu and Ehret (1990) is expected to result in 
significant savings of processing time and computer storage. 

In Chapter II, the WKB formulae are developed. The four types of normal modes 
are discussed. In Chapter HI, the normal modes and their eigenvalues for three abstract 
sound speed profiles are solved exactly. The WKB results are then compared to the exact 
solutions for an error analysis. Conclusions are given in Chapter IV. 
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II. WKB NORMAL MODES 



A. FIRST ORDER WKB APPROXIMATION 
1. The Mathematical Problem 

The WKB approximation to the solution of a dilTerential equation is an 
asymptotic expansion method applicable to a large range of equations. It is named after 
Wentzel, Kramers and Brillouin who used it separately but at about the same time in 
1926. However the principle of this technique was developed by Liouville and Green in 
1837. It was also used by Rayleigh (1912), Gans (1915) and Jeffreys (1924), among oth- 
ers. The WKB method is also known as the Liouville-Green or WKBJ approximation 
(the letter J is used to honour Jeffreys' contribution). 

To compute the acoustic normal modes the equation to be solved is the 
Helmholtz Equation 



dz 



( 1 ) 



where 



2 

2 CO 2 
= — ~^n 
c 

is the vertical wavenumber, Z„ is the «''' normal mode, z is the vertical coordinate, co the 
acoustic angular frequency, c = c(z) the sound speed and the horizontal wavenumber. 
With a pressure release surface and a rigid bottom the appropriate boundar>’ conditions 
are 



Z„ = 0 err the surface (2) 

dZ„ 

— ■ — = 0 at the bottom . (3) 

az 

Normal modes subjected to general boundary conditions are discussed in Appendix B. 
Note that, in (1) through (3), we have supressed the dependence in range which has no 
effect on the local modes. Equation (1), together with (2) and (3), is a Sturm- Liouville 
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problem where the Z„'s are the eigenfunctions and the k„’s are the eigenvalues. The 
eigenfunctions can be normalized using a normalization constant C„ such that 

(4) 



where 



I 



Z„Z„dz = <5„„ 



( 5 ) 



Note that the integration is over the entire depth in (5), 5^ is the Kronecker Delta, and 
Z„ and Z„ are orthogonal for n ^ m. 

There are several ways to obtain the first order WKB solution. A physical ap- 
proach is to consider the transmission of plane waves through a layered medium. The 
WKB solution is obtained as we let the the layer thickness approach zero. This physical 
approach will be discussed next. 

2. Physical Approach 

In each isospeed layer the solution to the Helmholtz Equation is given by 



Z„oce 






where k{„ is the vertical wavenumber in the layer. The transmission coeflicient from 
layer j to layer j+ 1 is given by (Kinsler et al., 1982) 









or 



^V+I 



1 



1 + 



1 

'^2n 



where = Kit'— k:^,. With 



^KL 



• 1 we get 



^7+1 = 



1 

^tn 



By letting | 1 = 1 we have 
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.y-1 



j-\ 



7^ = 



n = ^ 2 X! <r 

/ '"=• 






J 

= 1 



where 2> is the value at the interface of the /* layer and the (/ 4- I)'* layer and Az„ is the 
thickness of the m‘^ layer. Taking the thickness of the layers to zero we have 



Uz) = e 2J 



or 



\Xzn 



where is the vertical wavenumber at 2 = 2 , and 



0 



= Jj«„ 



\(iz . 



The vertical wavenumber generally can take on real or imaginary’ values. If 
K„ is real, i.e., > 0, the solution is 



Zn(z) = 



a’e‘‘^+b'e-‘^ 

\/Xzn 



or 



Zni^) 



a cos 4>+b sin 0 



(6) 



where a' and b' or a and b are constants to be determined by normalization and the 
boundary conditions. On the other hand if k„ is purely imaginary, i.e., < 0, we have 



Zniz) 



ae‘^+be 



( 7 ) 
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These last two formulae are precisely the WKB solution to the Helmholtz Equation in 
regions not close to a point where = 0 (i.e., a turning point). 

3. Basic Formulae 

By substituting the WKB solutions (6) or (7) in the Helmholtz Equation, one 
can easily show that they are exact solutions if 



Therefore, if |r( 2 )|<|i the W'KB solutions (6) and (7) are good approximations to the 
exact solution for k^„ > 0 and < 0, respectively. In Appendix A a more detailed and 
mathematical derivation of the WKB solution and its validity is presented. Let us call 
these solutions the first order WKB approximation. But (6) and (7) are not valid when 
= 0 (i. e., at a turning point) or even in a region where k^„ ~ 0. We need a solution 
valid at and near the turning point (i.e., in the critical region) to make the liaison be- 
tween the oscillatory region (kJ„|> 0) and the exponential region (k^„<^0). 

Let us suppose that there is a turning point at z = 0 and > 0 for z>0 and 
< 0 for z < 0. Consistent with a first order approximation, let us assume that near the 
turning point k]„ = yz where 





With a change of coordinate 




the Helmholtz Equation becomes 




which is the Airy Equation with the general solution given by 



Z„{0 = aAii0+bBii0 . 



The Airy Functions can be expressed as, with C > 0 , 



7 



BK-Q - 



AiiO = 




z?/(0 = 




/-„3(|ci-^„3(|c''0 

/-,;3(|<’'“)+'„3(|c’«) 



where J's are Bessel Functions and I's are Modified Bessel Functions. 

Now we have a complet set of first order approximate solutions covering the 
entire range of k^„. Summarizing, we have: 

(a) in the oscillator^' region (k^^O) 



Z„(z) 



a, sin ^+^1 cos 0 



\ 




( 8 ) 



(b) in the critical region (k]„ ~ 0) 



X:n = V- 



r 1/3 

!^ = -y ‘ Z 



Z,{z) = a^Ai{!:)+b^Bm , 



(9) 



(c) in the exponential region (k^„<^0) 



2„(2) 



a^e 

J I ><7n 1 



( 10 ) 



To assure continuity in Z„(z) at the boundaries between regions, we must relate 
the constants a, and b^ carefully. An asymptotic expansion of the critical region solution 
for large positive values of k]„ (or — C > 1) must match the solution in the oscillatory re- 
gion and for large negative values of (or C > 1) niust match the solution in the expo- 
nential region. The connection formulae between the three regions are given by 
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■ “ 7 = sm{4>+ -J ) 

V • I \ ^zn 

—zA=re* *-oBi{Q -* J — cos(0+-j) 

\' I ^2n I n/ ^2n 

where the left sides match the Airy Functions for k]„-4.0 and the right sides match for 
(Bender and Orszag, 1978) with 

r~ - 1/6 

(7 = ^^/ 7T y 

Therefore, (8), (9) and (10) can be recast respectively as 

(a) oscillatory region solution 



a sin(0+ )-\-b cos{4>+ ) 



(b) critical region solution 

Z,{z) = aaAtX0+b<^Bi(0 , 



(c) exponential region solution 



Zni^ 



(a/2)c ’^+be'^ 

\i I I 



(ll.cz) 



(ll./z) 



(ll.c) 



The constants zc„, a and b are determined by normalization and the boundary 
conditions. The equations for the eigenvalues k„, i.e., the characteristic equations, in- 
volve the solution in the oscillatory region (11. a), as will be discussed later. 

Equations (ll.a,b,c) are not easy to use in the computation of normal mode 
shapes. Where does one stop with one formula and start with another? This problem is 
avoided if, after the determination of the eigenvalue and the constants a and b, the 
computation of the normal mode is done using the Langer Formula (Nayfeh, 1973; 
Bender and Orszag, 1978) 



Z,(z) = 2^/7(-|-|0|) 



>/ 6 . 



aAi 






( 12 ) 
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with 



and 



(J) = 




K,„ I dz 



s = -zl\z\ . 

This formula has the advantage of giving a single and continuous solution for the entire 
range of k^. Bender and Orszag (197S) have shown that for ~ 0 and K^„ = yz the 
formula gives exactly {1 l.b), for k^I> 0 it asymptotically approaches (11. a) and for 
it approaches (1 l.c). 

In rectrospect, in order to obtain the first order WKB solution, an algorithm 
must include the following steps: 

(a) application of a boundary condition to (ll.a,b,c) to get the relationship 
between the constants a and b\ 

(b) application of the other boundaiy condition to (1 l.a) to get the character- 
istic equation; 

(c) with (12) compute Z„(z); 

(d) with (4) and (5) normalize Z„(r). 

4. Comments 

The WKB solution is a high frequency approximation. This means that it gets 
better as the frequency gets higher. It must be noted that, although Normal Modes are 
a full-wave exact solution to the wave equation, the WKB modes are approximate sol- 
utions and their accuracy is frequency dependent. 

B. DETERMINATION OF WKB MODE PARAMETERS 

Now that we have a uniformly- valid first order solution to the Helmholtz Equation, 
we are ready to apply the boundary conditions to get expressions for the k„'s and the 
constants a and b (or their equivalents). There are four classes of normal modes to 
consider. Each class has different mathematical expressions for the mode parameters (a, 
b , or their equivalents, and k„). Therefore, for an arbitrary sound speed profile the first 
procedure for normal mode calculation is to determine which class each mode falls into 
and then go through the steps described at the end of the previous section. 
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1. General WKB Forimilae 

We start by deriving the more general formulae which are applicable to modes 
that do not have turning points in close vicinity of the boundaries. 
a. Class I: Pressure Release - Turning Point 

First let us consider a mode whose oscillatory region is bounded by a turn- 
ing point at the depth z = z and the surface (z = 0). With the bottom at z = //, the 
boundary conditions are 



Z„(0) = 0 



(13) 



dz 



{H) = a . 



(14) 



We will call this class of modes as PR-TP (Pressure Release - Turning Point). 

In the exponential region (z<|z < H), the WKB solution can be recast as 



Z„ oc 



cosh(0-0fc) 

V I Xzn I 



(15) 



with 



4> = 




In order to satisfy the boundary condition (14), 4>i, must be given by 



<!>,= 




I dz— tanh 





dz 



We need to connect (15) to the solutions in the critical and oscillatory regions next. After 
application of the connection formulae, the results in the three regions are: 

(a) exponential region 




cosh(<^-0^) , 



(IM 



(b) critical region 

= v(2 -z) 
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C = -y'^'(£-z) 

Z„ = eKAi{0+ oBiiO , 

(c) oscillatory region 




{I6.b) 



(16.C) 



with 



4 > = 




K,„ I dz 



in this region. The corresponding Langer Formula that asymptotically matches 
(16.a,b,c) is 




with 



j = 



A 



z— 2 



Z-Z 



and 



0=1 I I dz. 

As the oscillatory region solution must satisfy the surface boundary condi- 
tion, we have 






cos 





sin 




= 0 



(18) 



or 
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^ ( 1 ^ 

dz=[n — - In:— tan ( — — 

V 4 y 

This is the characteristic equation for PR - TP normal modes. The solutions to the 
characteristic equation give the eigenvalues k„'s. 

b. Class II: Turning Point - Rigid Bottom 

The next class of normal modes to be considered has the oscLllatorv' region 
between a turning point and the bottom. It will be called TP-RB (Turning Point - Rigid 
Bottom). For this class, let us express the solution in the exponential region, 0 < z<^z, 
as 





Z„ oc — , sinh((i>— d>r) 



with 



d> = 




dz 



and 



4>s = 




Using the connection formulae, we get for the three regions: 
(a) E.xponential Region 

?===- smh{4>-(f),) 

yj I I 



(b) Critical Region 



4 = >’M) 

C = -y'/'(2-z) 
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(c) Oscillatoiy Region 






= -= sin('*+ T )- "<■<'+ f ) 

V ^zn ^zn 



(19j 



with 



4>=j\K^„\ dz 



in this region {z4.z < H). The appropriate Langer Formula for this class is 




and 






\dz. 



Equation (19) must satisfy the bottom boundary condition (14), therefore, we have 



where 



sin 



I 4 



+ tan 






( 21 ) 



Z) = 



I 

Ik 



U=H 



( 22 ) 



It follows from (21) that the characteristic equation for TP - RB is 
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K— tan 



e‘^’+ 



£ 

2 



e 












-De*‘ 



for D <0, and 



K:ndz = 




for D>0. 

c. Class III: Turning Point - Turning Point 

The third class of modes has the oscillatory region between two turning 
points, at z = z, and z = Zj with z, < Zj. They will be refered as TP-TP (Turning Point - 
Turning Point). 

Let us express the solution in the exponential region near the surface, 

0 ^ z<Zi, as 



2 . = - 



V I I 



sinh((ji), -4>,) 



where 



and 



A 





and the solution in the exponential region near the bottom, Z 2 <^z < H, as 



Z„ = cosh(0^-(/)*) 

V I Kr J 



with 
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and 



4 > 



^ _ 
2 




Xin I 




In addition, let us express the solution in the critical region centered at z = z, as 

Z„ = /!/(:,)- -y (23) 

with 

2 / ^ \ 



and 

Cl = -yj^^(z-zi), 

and in the other critical region, centered at z = Zj, as 

Z„ = C2C^‘(T2/1/(C2)+ ^ (24) 

with 

xl„ = y2i^x-^) . 

<72 = sjn 

and 

Cz = -y’2\^2-2)- 

By letting 
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Cl = 






- 24 >, 



and connecting (23) to the solution in the oscillatory region i,<^z<|z 2 . get 



where 



2„ = - 



V ^zn 



sin((/)t+-j-a,) 



( 25 ) 



4> 



r=jVzj 



dz 



and 



ot) = tan 



-1/ c 



-‘t’. 






Note that, in obtaining (25), the trigonometric identity 

a cos 6+b sin 6 = sja'+b'^ sin^0+ tan”' ^ 
has been used. In the same way, by letting 



(26) 



C') — ' 






-2<Pb 



and connecting (24) to the oscillatory region solution we get 

Z„ = — J=- sin(0^+ -j +«2) 

where 



(27) 



4> 



2 = f Vrn I 



dz 
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and 




Realizing that the two solutions, (25) and (27), in the oscillatory region 
must be identical, we obtain 

sin(</)f+^-a,) = sin(07+-^+a2) . (28) 

Using (28), the characteristic equation for this class of normal modes is obtained: 

rt— ^TT+a,— «2 . 

To compute the normal mode w'e can use (20) for z and (17) for 2 > z^. 
d. Class IV: Pressure Release - Rigid Bottom 

Finally, we must consider modes having no turning points. We call this class 
PR-RB (Pressure Release - Rigid Bottom). Here we express the solution, which is always 
oscillator>‘, as 




Z„oc 



sin 4>-\ cos <f) 

Jxzn 



(29) 



with 



<^ = [ \X;n\d2 ■ 

•’0 

Since the surface boundary condition must be satisfied, we must have b = 0, and hence 



2 .= 




sin 0 . 



(30) 



On the other hand, (30) must satisfy the bottom boundary condition (14) and tliis leads 
to the characteristic equation equation for the PR - RB modes: 
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rw 

K^„dz = nn— tan 

Jq 

for D>0, and 



rH 

K^ndz = rt7T+ tan 

•*0 

for Z) < 0, where D is given in (22). 

2. Formulae Associated With Near-Boundary Turning Points 

In our derivation of the formulae in the previous section, we have applied the 
boundary conditions to the oscillatory or exponential region solutions, assuming that the 
boundaries are nowhere close to any turning point. In the special case that a boundary 
lies inside a critical region, the respective boundary condition must be applied to the 
critical region solution instead. Different formulae associated with the first three classes 
of modes for this special case will be derived next. 
a. Class I 

PR - TP modes have lower turning points. Since the solution in the critical 
region must now satisfy the bottom boundaiy condition we obtain 





where 



and 




Ai' and Bi' are the derivatives of the Airy Functions with respect to C and are defined 
by (with C ^ 0) 



,di'(-0 = -yC 
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V3 






/i/'(0 = -yC 



Bi'iO — 



i-,„{ie'^)-hi,{j-i’'^) 



Defining 



A, = 5/'(C;,) 



and 



B, = Ai'i^f,) 



and applying the connection formulae, the solution in the oscillatory region becomes 



\'^zn 



5. 



V '^zn 



An application of the surface boundary condition gives the follouing characteristic 
equation: 






b. Class II 

TP - RB modes have upper turning points. To satisfy the surface boundary 
condition, we require the solution in the critical region to vanish at the surface, i. e., 

Z„ = Bi{i:s)oAi{0-Ai{i:s)oB!{0 



with 



Cs = y 



1/3a 



and 
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Defining 



y = 



dz 




= Biyls) 



and 



— Ai(Cs) 



and applying the connection formulae, the solution in the oscillatory region becomes 




The subsequent application of the bottom boundary condition leads to the following 
characteristic equation ; 




for D <0, and 




for D >0, where D is defined in (22). 
c. C/ass III 

With the above results, we can easily that only two modifications in the 
general formulae for TP - TP modes are required. These include replacing e‘^> by and 

by B 2 , if the upper turning point is close to the surface, and by Ai and — 

by Bi, if the lower turning point is near the bottom. 

3. Criterion For Formulae Selection 

We must define a criterion for switching from the general formulae to the for- 
mulae associated with near-boundary turning points. It was found by Bender and 
Orszag (1978) that the critical region extends on each side of the turning point until 
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I Cl ~ 1 . In accordance, we will use the general formulae unless the distance between a 
turning point and a boundar>' corresponds to values of ICl smaller than 1. 

4. Normalization 

Once normal mode parameters are computed using the appropriate formulae 
including the characteristic equations, the normalization of the normal modes can be 
achieved by numerical integration over depth. The normalized modes (Z„'s) are related 
to the Z„ s by 

Z„ — c„z„ 

n n^n 



with 




('ll 

Zldz . 

*'0 



5. Parameter Storage 

We now define the necessar\- parameters required to completely characterize 
each mode. The first parameter is obviously the class number. By checking the formulae 
developed in the previous sections it is seen that all the constants (0„ 0*, D, etc.) can 
be computed from the horizontal wavenumber, the depths of the turning points and the 
depth of the ocean. Therefore, the parameters required to parameterize each class of 
modes are: 

a. Classes I And II 

1. Class number 

2. Horizontal wavenumber 

3. Normalization constant 

4. Depth of the turning point 

b. Class III 

1. Class number 

2. Horizontal wavenumber 

3. Normalization constant 

4. Depth of the upper turning point 

5. Depth of the lower turning point 
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c. Class IV 



1. Class number 

2. Horizontal wavenumber 

3. Normalization constant 
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III. ACCURACY TEST 



A. EXACT ANALYTICAL SOLUTIONS 

To quantify the accuracy of the first order WKB approximation, we will compare 
the WKB solutions to exact analytical solutions to the Helmholtz Equation for three 
abstract sound speed profiles. The exact solutions are derived in this section. We use 
as boundary conditions pressure release surface and rigid bottom for all three cases. 

1. Positive Exponential Profile 

In this upward refractive profile, sound speed increases exponentially with the 
depth and is given by 



where /? is a constant. The surface is at z = 0 and the bottom at z = H. The Helmholtz 
Equation is 



c(z) = 




or 





( 31 ) 



with kI = . 



With a change of variable 




(31) becomes 




After division by and using 




( 32 ) 
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and 



a 



n 






(33) 



we obtain 



dZ„ 



2 






This is the Bessel Equation and its general solution is 

Zn = (^JA«oX)+bY M qx) 



The boundar>’ conditions are 



Z „(2 = 0)^Z„(;r=l) = 0 



dZ^ dZf^ 



The application of (35) and (36) to (34) results in two algebraic equations 

(®0Xh) = 0 . 



(34) 

(35) 

(36) 

(37) 

(38) 



As this linear system of algebraic equations is homogeneous, there is a nontrivial sol- 
ution only if the corresponding Jacobian is zero, i.e.. 

This characteristic equation must be solved in order to obtain the a„'s, which are the 
scaled eigenvalues. After the a„'s are evaluated, the k„'s can be determined using (33). It 
follows from the surface boundary condition, expressed in (37), that the specific solution, 
aside from a multiplicative constant, is 

z,(2),>= ■ 

The constant of proportionality is given by normalization. 
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2. Negative Exponential Profile 

The second profile considered is a downward refractive one. Here, sound speed 
decreases exponentially with depth and is given by 

c{z) = c^e~^' • 



Thus, the Helmholtz Equation is 

= 0 (39) 

dz 

1 CO 

where Kj = — . 

Co 

With 



<? = 



e 






and tto and a„ given by (32) and (33) respectively, (39) can be recast as 



LA. 




(40) 



The general solution to (40) is 



2„ = n./Jao(f)+*FJaoO • 



Using the same anology as in the previous case, w'e obtain as characteristic 
equation 

- 0 (41) 

and specific solution 

Z,(z)c^ • 



The K„'s can be found using (33) after solving (41) for the a„'s. 

3. Hyperbolic Profile 

In this third case we assume that the sound speed has a minimum at 2 = Zj. The 
sound speed profile is 



c(z') = Cq cosh(^z') 
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with 



z = z—i 



The Helmlioltz Equation to be solved is now 






dz 



,2 






(t) 



2 1 






cosh^(^z') 



Z, = 0 . 



(42) 



With the following change of coordinates 

X = tanh(^z') 



(42) becomes 



/?'(l-x')-^ 



i\-x) 



2v 



dx 






(43) 



After dividing (43) by ^^(1-y^) and using the definitions given in (32) and (33), (43) can 
be recast as 



dx 



n 2,^^. 



+[ an- 



2 ^-n 



0 2 
1-x ] 



z„ = 0 . 



(44) 



This is the Associated Legendre Equation. The general solution to (44) is 

Z, = an"(x)+^a"(x) , 



with 



and 

v(v+l) = ao . 

The functions Pt{x) and 2f(x) are the Associated Legendre Functions of the first 
and second kind of order n and degree v . Applying the boundary conditions we get the 
system of equations 

<i'’;'(x5)+i’2.''(;!j)-o (J5) 
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aK{y.u)+bQ":{xu) = ^ 



with 



X5 = tanh(-^2o) 



and 



;>;// = ianh[^(H-zo)] • 

Thus, the corresponding characteristic equation is 

= 0 . (46) 

Solving (46) for m gives the eigenvalues a, 's and hence the horizontal wavenumbers 
K„ 's. Using (45) the solution becomes 

Z„{i') Q:\Xs)P:\ tanh m-P:\xs)Q:\ tanh pz') . 

B. ANAL\TICAL EVALUATION OF WKB PHASE INTEGRALS AND 
DERIVATIVES 

Now we use the WKB formulae to obtain the first order WKB solution for the three 
analytical sound speed profiles. The objective in this section is to evaluate analytically 
all the related WKB integrals and derivatives so that numerical errors can be eliminated 
in one of the error analysis. Horizontal wavenumbers computed using both the exact and 
approximate (i.e., WKB) characteristic equations will be compared in the next section. 
The evaluation of the horizontal wavenumbers is emphasized because they are the key 
parameters in normal mode computations. 

1. Positive Exponential Profile 

The vertical coordinate z is positive downward, the surface is at z = 0 , the bot- 
tom at z = // and the sound speed profile is 

C(z) = . 



With 



'^0 = 7 ^ 

the corresponding wavenumber profile is 
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K — Kq€ 

Let us first consider modes that have a turning point at z = f . for 0<z<z 
(oscillatory region), is given by 

/ -2/37 T 
= ^ -Yn 



with 



Kq 

The spatial phase (f> in this region is given by 



0 = 




K,„\dz . 



(47) 



With the chanse of coordinate 



(47) becomes 



-2(32 

X = e 






y„ tan 



-1 



y-Vn 



nX(^) 






-sJX-Yn 



■^x(^) 



In the region z -^z ^ II (exponential region), we must use 



/ 2 2-2pz 



I I III 

I ^<2n I = V 



The spatial phase is now 



= I J ^zn I dz , 



or, with 
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We must also compute 



D = 



d\K, 



2 1 K , J 



2 dz 



h=H 



which can be equated as 



D = 






-HH 






Figure 1 shows a PR - TP mode in this upward refractive medium. 
If the normal mode is PR-RB the phase is given by 



<^» = I \^2n\dz 
♦’o 



or 




Vn tan 




X{2) 

X(0) 



with 



y.{z) = e 



and 



n 4l> 

For each case we must also apply the respective characteristic equation. 

2. Negative Exponential ProFile 

For this profile we have two types of normal modes: TP - RB and PR - RB. 
The sound speed is 
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Figure 1. PR - TP Mode 



c{z) = CQe . 

Consequently, for a TP - RB normal mode and for 2 > z (oscillatory' region) we have 

= KqJ . 



So in this region the phase is 
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«0 






\ A fn 



/ n 



tan 










with 



2p2 

X = e ^ 



Also we have 



kIH 

2 



In the exponential region (0 < z < z), 



\>^2n\= «0\/ y'n 



—e 



Ijiz 



and thus the phase is 





A,2_<c2 

V 'n - 










(48) 



with 



^ = . 

Figure 2 shows a TP - RB mode in such a downward refractive medium. 
For a PR - RB normal mode we have 



(j) = 







-Vn tan 




-i;tW 

-*x(0) 



and D is given by (48). Again, for each of the classes the respective characteristic 
equation must be used. 

3. Hyperbolic Profile 

We consider two types of normal modes: TP - TP and PR - RB. For the type 
TP - TP and in the oscillatory' region (z', < z < z'j) the phase will be 
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(f) = 




( 49 ) 



Using 



^ = tanh(^z') , 
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and 



(49) becomes 



^0 

P 



2 _ I JfiL 

Hn ^ 2 

»<0 



0 = «o 



j — ^ 

• -i/ ^ ^ • ->/ M 

sin (^1^ j-sjl-rin Sin 









For the exponential region near the bottom {H>z'> z'j) the phase is 



4 > = 




^ zn J 



or 



«n 


lr> 




1 « 1 1 


\+t yn\! ^^-nl-{^-V)+yl 


-l 


2 


in 


^ + \J 


> rn 


^ yn\j ~(l+'^)+)'n 


- 



where 




For this type of normal mode D is given by 

■ 2[4-4(i-<i,)F ' 

In the exponential region near the surface (z', > z') the phase is 
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A 



If the mode is PR - RB the phase will be 










and D is given by 









Figure 3 shows a TP - TP mode trapped by the refractive index. 

C. COMPARISON OF RESULTS 

1. Horizontal Wavenumber Error Analysis 

The exact characteristic equations derived in section A for the three abstract 
profiles were solved using iterative procedures to obtain the benchmark k„ 's. Similarly, 
the approximate WKB characteristic equations developed in the previous chapter were 
also solved for the three profiles with the use of the algebraic equations developed above. 
Very low frequencies were chosen for the comparison intentionally. This would give the 
WKB method a real test, since WKB is inherently a high-frequency approximation. For 
the exponential profiles, we used &>=10 or a frequency of 1.59 Hz. For the 
hyperbolic profile the frequency used was even lower, it was 0.23 Hz. Some of the 
computed horizontal wavenumbers (k„) are presented in Tables 1 through 5. The unit for 
K„ is inverse meter (w‘). The depth of the ocean is taken to be 10000 m. 

In Tables 1 and 2 the exact and WKB results, as well as the absolute and relative 
errors, for the positive exponential profile are displayed. Tables 3 and 4 show the results 
for the negative exponential profile. In Table 5, the exact and appro.ximate results for 
the hyperbolic profile are compared. At 0.23 Hz, mode 2 is TP - TP and mode 3 is PR 



Overall, we can see that the absolute error (the absolute value of the WKB re- 
sult minus the exact result) varies between 10"* and lO*’ nr' and the relative error (the 



- RB. 
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Figure 3. TP - TP Mode 



absolute error divided by the exact value) varies between 10‘* and 10'^ It is noted that 
for the modes with turning points not in close vicinity to the boundaries the absolute 
error usually stays below 10'^ »r'. The fact that the error increases when a turning 
point is close to a boundary is easily explained. The solutions in the critical regions are 
the less accurate because they use a linear approximation for So when they are used 
to match the boundary conditions the error increases. 
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Table 1. POSITIVE EXPONENTIAL PROFILE, PR - TP MODES 



.Mode 

# 


Exact K„ (/«-') 


WKB (m-‘) 


Absolute 
Error (tn-') 


Relative 

Error 


5 


.43861 X 10-^ 


.43866 X 10-2 


5 X 10-’ 


1.1 X lO "' 


6 


.40441 X 10-2 


.40445 X 10-2 


4 X 10-2 


1.0 X lO-* 


7 


.37227 X 10-2 


.37230 X 10-2 


3 X 10-2 


.8 X lO--* 


8 


.34179 X 10-2 


34182 X 10-2 


3 X 10-2 


.9 X 10-^ 


9 


.31274 X 10-2 


.31277 X 10-2 


3 X 10-2 


1 X 10-“ 


10 


.28556 X 10-2 


.28568 X 10-2 


1.2 X 10-‘ 


4.2 X 10-“ 



Table 2. POSITIVE EXPONENTIAL PROFILE, PR - RB MODES 



.Mode 

M 


Exact K„ (m-‘) 


WKB K„ (A?r‘) 


Absolute 
Error (m"') 


Relative 

Error 


12 


.23190 X 10-2 


.23041 X 10-2 


1.5 X 10-5 


6.3 X 10-2 


13 


.18524 X 10-2 


.18490 X 10-2 


3.4 X 10-' 


1.8 X 10-2 


14 


.10756 X 10-2 


.10736 X 10-2 


2 X 10-' 


1.9 X 10-2 



Table 3. NEGATIVE EXPONENTIAL PROFILE, TP - RB MODES 



.Mode 

rr 


Exact K„ (/?r‘) 


WKB K„ (/n-‘) 


Absolute 
Error (m-') 


Relative 

Error 


23 


.82307 X 10-2 


.82308 X 10-2 


1 X 10-’ 


1 X 10-2 


24 


.79463 X 10-2 


.79464 X 10-2 


1 X 10-2 


1 X 10-2 


25 


.76664 X 10-2 


.76665 X 10-2 


1 X 10-2 


1 X 10-2 


26 


.73903 X 10-2 


.73904 X 10-2 


1 X 10-2 


1 X 10-2 


27 


.71158 X 10-2 


.71155 X 10-2 


3 X 10-2 


4 X 10-2 



In general, for each type of mode the WKB approximation gets better as mode 
number increases. The only exception is when turning points are close to the boundaries. 
To see this, let us examine the results for the positive exponential profile (Tables 1 and 
2). Starting from the lowest modes that have one turning point, the error decreases as 
mode number increases. At mode number 10, the accuracy of the WKB method de- 
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Table 4. NEGATIVE EXPONENTIAL PROFILE, PR - RB MODES 



Mode 

r? 


Exact K„ {nr'^) 


WKB K„ {nr') 


.Absolute 
Error (am*') 


Relative 

Error 


29 


.65250 X 10-' 


.65091 X 10-' 


1.6 X 10-' 


2.4 X 10-' 


30 


.61746 X 10-' 


.61683 X 10-' 


1.6 X 10-' 


1.0 X 10-' 


31 


.57723 X 10-' 


.57688 X 10-' 


3.5 X 10-' 


6. 1 X lO--* 


32 


.53082 X 10-' 


.53058 X 10-' 


2.4 X 10-' 


4.5 X lO"' 


33 


.47670 X 10-' 


.47652 X 10-' 


1.8 X 10-' 


3.8 X 10^ 



Table 5. HYPERBOLIC PROFILE 



M ode 


Exact (nr') 


WKB (am-') 


Absolute 
Error (m-') 


Relative 

Error 


2 


.901 14 X 10-' 


.89881 X 10-' 


2.3 X 10-' 


2.2 X 10-' 


3 


.74139 X 10-' 


.72953 X 10-' 


1.2 X 10-' 


1.6 X 10-' 



creases because the turning point comes too close to the bottom boundary. For mode 
12 and higher a turning point does not exist and the oscillatory region is bounded by the 
physical boundaries. After the change of mode class the error starts to decrease again 
as mode number increases. 

Errors in the k„'s for all cases are small enough that they do not significantly 
influence the vertically integrated phases 0 because the maximum ocean depth is of the 
order of km's. In the horizontal direction the horizontal phase is approximately K„r (this 
phase is exact when c = c(z)). The error in k„ limits the range for which the use of WKB 
is accurate. If we want to keep the phase error below a few degrees, the tolerable error 
in K„ is of the order of 10”^ w* for a range of 10 km, 10*‘ w' for 100 km and 10"^ m*' 
for 1000 km. All the k„'s associated with the modes that do not have interactions with 
the bottom boundaiy, as calculated from the WKB method, have errors less than 
10"‘ nt~\ implying that the method is good to at least 100 km at 1 Hz. For higher fre- 
quencies the results would be much better. 

2. Errors In The Interference Distances 

The acoustic pressure square at far field can be written as (assuming no range 
dependence) 
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oo 



o)Z„(--)+ 



;i=l 



oo 

NT' 

/ COS[(K^-K„)r] 

«=1 

m^n 



where the A„'s are constants and the acoustic source is at z = Zq and r = 0 (Clay and 
Medwin, 1977). We can see that the interference terms are functions of the differences 
of horizontal wavenumbers. The interference distance (or interference wavelength) de- 
fined by 



A 



nm 



2n 

Km-^n 



is therefore an important parameter for transmission loss calculation (Chiu and Ehret, 
1990). In general, dominant interferences are between adjacent modes (Chiu and Ehret, 
1990). So let us see what is the size of the errors in Ak„ = kJ computed using the 
WKB approximation. In Tables 6 and 7 we display the exact and WKB Ak,'s, as well 
as the absolute errors. As expected, we can see that the errors in Ak„ var\' in a similar 
way as in k„ and are slightly smaller. 



Table 6. POSITIVE EXPONENTIAL PROFILE, PR - TP MODES 



n 


Exact Ak„ (att') 


WKB ^K„ (m-') 


Absolute Error 
(/«-') 


5 


.3420 X 10-^ 


.3421 X 10-3 


1 X 10-’ 


6 


.3214 X 10-3 


.3215 X 10-3 


I X 10-’ 


7 


.3048 X 10-3 


.3048 X 10-3 


0 


8 


.2905 X 10-3 


.2905 X 10-3 


^ 0 


9 


.2718 X 10-3 


.2709 X 10-3 


9 X 10-’ 



3. General Numerical Method 

The WKB results in the previous section were obtained using analytical means 
rather than numerical methods. The reason for that was we wanted to see how accurate 
the WKB approximation is regardless of the numerical methods used to evaluate inte- 
grals and derivatives. As we see, absolute errors in the order of 10"'’ m*' for the k„'s are 
achieved at 1 Hz. This means that the WKB solution is very accurate, especially because 
we used very low frequencies. 
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Table 7, POSITIVE EXPONENTIAL PROFILE, PR - RB MODES 



n 


Exact Ak„ (nr‘) 


WKB Ak„ (m->) 


Absolute Error 
(m->) 


12 


.4666 X 10-' 


.4551 X 10-' 


1.15 X 10-' 


13 


.7768 X 10-' 


.7754 X 10-' 


1.4 X 10-' 



For an arbitrary profile all the integrations and difTerentiations must be done 
numerically. It is expected that the errors will increase due to numerical noise. The 
Fortran program WKBGEN (Appendix C) developed in this thesis can be applied to an 
arbitrary profile. For each mode it finds class, horizontal wavenumber and depths of the 
turning points. This program was also applied to the three abstract profiles. Some nu- 
merical results for the k„'s are presented in Tables 8 through 1 1. 

Tables 8 and 9 show' the exact and numerical WKB results and the respective 
absolute errors for the positive exponential profile. Tables 10 and 11 show to the cor- 
responding results for the negative exponential profile. 



Table 8. POSITIVE EXPONENTIAL PROFILE, PR - TP MODES 



Mode 

u 


Exact 


WKB K„ (ni->) 


Absolute 
Error (m-‘) 


5 


.43861 X 10-' 


.43874 X 10-' 


1.3 X 10-' 


6 


.40441 X 10-' 


.40454 X 10-' 


1.3 X 10-' 


7 


.31221X 10-' 


.37234 X 10-' 


7 X 10-" 


8 


.34179 X 10-' 


.34184 X 10-' 


5 X 10-’ 


9 


.31274 X 10-' 


.31284 X 10-' 


1 X 10-' 


10 


.28556 X 10-' 


.28574 X 10-' 


1.8 X 10-' 



As w'e can see, the errors in the k„'s computed using the numerical method are slightly 
larger than the previous results (Tables 1 through 4) but are approximately in the same 
order of magnitude. 
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Table 9, POSITIVE EXPONENTIAL PROFILE, PR - RB MODES 



Mode 

ii 


Exact K„ (m*‘) 


WKB (a«-') 


Absolute 
Error (w-‘) 


12 


.23190 X 10-' 


.23044 X 10-' 


1.5 X 10-5 


13 


.18524 X 10-' 


.18494 X 10-' 


3.0 X 10-« 


14 


.10756 X 10-' 


.10744 X 10-' 


1.2 X 10-‘ 



Table 10. NEGATIVE EXPONENTIAL PROFILE, TP - RB MODES 



Mode 

a 


Exact K„ (m-‘) 


WKB (m->) 


Absolute 
Error (/u-') 


23 


.82307 X 10-' 


.82309 X 10-' 


2 X 10-' 


24 


.79463 X 10-' 


.79469 X 10-' 


6 X 10-’ 


25 


.76664 X 10-' 


.76669 X 10-' 


5 X 10-' 


26 


.73903 X 10-' 


.73909 X 10-' 


6 X 10-’ 


27 


.71158 X 10-' 


.71159 X 10-' 


1 X 10-' 



Table 1 1. NEGATIVE EXPONENTIAL PROFILE, PR - RB MODES 



Mode 

# 


Exact K„ (au-') 


WKB K„ (m-') 


Absolute 
Error (AArO 


30 


.61746 X 10-' 


.61689 X 10-' 


5.7 X 10-5 


31 


.57723 X 10-' 


.57689 X 10-' 


3.4 X 10-5 


32 


.53082 X 10-' 


.53059 X 10-' 


2.3 X 10-5 


33 


.47670 X 10-' 


.47649 X 10-' 


2.1 X 10-5 



41 



IV. CONCLUSIONS AND RECOMMENDATIONS 



The WKB approximations of acoustic normal modes seem to give results accurate 
enough to be of practical use. Although WKB is inherently a high frequency approxi- 
mation, meaning that it works better for higher frequency, the results from our tests 
show that, even for frequencies around one Hertz, this technique has an appreciable 
accuracy. 

When applying the WKB algorithm for arbitrary sound speed profiles the determi- 
nation of the class of each normal mode must be done carefully. Each class has different 
formulae. There are a total of four classes in a single duct or channel environment. 

A small weakness of the method is that the error in the WKB solution increases 
when a turning point is very close to a boundary. However, the corresponding errors are 
still very small for the exponential and hyperbolic profiles used in this study. 

The exact analytical solutions to the Helmltoltz Equation can also be used for 
comparison to any other methods. These exact solutions are subject to pressure release 
surface and rigid bottom boundary conditions. Exact solutions can also be found for 
any boundary conditions. Specifically, we can maintain the pressure release surface 
boundary condition, which is a good assumption, and use a non-rigid bottom boundary 
condition. 

Some difficulties were found when working with Bessel and Associated Legendre 
Functions. The IMSL subroutines used for Bessel Function evaluation cannot handle 
some orders and values of the argument. With respect to the evaluation of Associated 
Legendre Functions, subroutines are not available in the IMSL libraries. Some Fortran 
programs were coded to compute them. These programs also work only for certain 
ranges of order, degree and arguments of the functions. Further programming work 
concerning these transcendental functions is recommended. 

WKB formulae for modes trapped in the water column over a non-rigid bottom were 
derived in Appendix B although they are not implemented for this analysis. A general 
mathematical e.xpression for boundary conditions for normal modes is also presented in 
Appendix B. The use of this mathematical expression is not a trivial matter because the 
expression depends on the density and the sound speed profile of the sediment. Further 
studies on using WKB algorithms for arbitrary bottom boundary conditions are recom- 
mended. 
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In this thesis only single duct,'channel environments have been considered. Essen- 
tially, we have ignored double duct/ channel problems. However, the WKB approxi- 
mation method has been applied to other nonacoustic but equivalent double 
duct, channel problems (for example, transmission of electromagnetic waves through 
potential barriers) with sucess (Bender and Orszag, 1978). 

In conclusion, the WKB method allows for accurate and fast computation of normal 
modes and their eigenvalues. In addition, it allows for storage of the results in terms of 
only a few parameters for each mode. 
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APPENDIX A. FIRST ORDER WKB THEORY 



With 



d<t> = I dz 



the Helmholtz Equation becomes 



d^Z„ 1 d\K,J 

d(t>^ \k,„\ d<P 



dip 



+Z„ = 0 . 



(^. 1 ) 



This equation has some similarities with a Bessel Equation. Let us try a solution in the 
form 



Z„ oc 




F{P). 



Substituting this trial solution form in (A.l) we get 



d^F 1 dF 

dip^ <i> dp 



l+r( 2 )- 



( 1 / 2 )^ 

0' 



f =0 



(A.2) 



with 



-■w-i+l 




d\K,n\ 

dz 



1 1 d^\K,„\ 

2 dz^ 



If r(c) is negligible, i.e., 1 r{z) 1 <|1, (A.2) can be approximated as 



d^F \dF 


\ (1/2)' 


dp^ P dp^ 





(^. 3 ) 



which is a Bessel Equation. 

The general solution to (A. 3), for > 0, is 

FiP) = a'J^i2{P)+b'Y,i,iP) {A.4) 

where 7 , ,2 and l ’,;2 are the order 1/2 Bessel Functions of the F' and 2"‘‘ kind, respectively. 
But since 



44 



sin X 



' 1/2 



OC ■ 



v'X 



y 



1/2 



OC 



COS jr 



(A. 4) becomes 



F{<f>) = a 



sin 0 



+b 



cos 0 



and thus 



a sin 4>+h cos 0 

2„(2) = = 



If < 0, instead of (A. 3) we would get the Modified Bessel Equation whose general 
solution is given by 



/■(0) = a'/f,p(0)+/l'/,/2(0) 



where A ^,/2 and A/2 are the Modified Bessel Functions. Since 



A'i/2(0)oc-^ 
V ^ 



cx: — =r 

n/0 




we now have 



f(0) = 



ae ^+be'^ 



and thus 



Zn{2) = 



ae_^_+b^ 

V U7 J 



45 



Let us assume that there is a turning point at z = 0, k^„>0 for z > 0, < 0 for 

2 < 0 and near the turning point k^„ = yz. So, close to the turning point we have 



4> = j |k,„Uz = 5y y'^^lzl 



3/2 



with s = z/ 1 z I , and after some algebra, 

1 d\K,„\ 1 



Using this in (A. I) we get 






1 dZ„ 

+ -:r, rf- +Z„ = 0 



d<t)^ 30 J0 

Let us now look for a solution near the turning point in the form 

z„oz(^^4>yih\<j>). 



Substituting this solution fonn in (A. 5) we get 



d^F 1 dF 

d<t>^ 4> d4> 






1 - 



(i/3r 



f =0 



{A.S) 



whose general solution is 



where B represents Bessel Functions. For k\„ > 0 they are y.j/j and J^y For < 0 they 
are /.,/3 and I^|y Relating Z„ with F we obtain 

(y 0)'''[qfi-,/3(0)+C25,/3(0)] 



or 



2^^ = V y U I 



C, 5-i/3( ± y U I +c2^i/3( ± y u I 



(,4.6) 
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where the plus sign is for k ]^ < 0 and the minus sign for k \^ > 0. A w^ay to avoid the in- 
convenience of sign switching is to use Airy Functions which are related to (see 
Chapter 111, Section A). 
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APPENDIX B. GENERAL BOUNDARY CONDITIONS 



A. GENERAL FORMULA 

Until now we have assumed that the boundary conditions are pressure release at the 
surface and rigid bottom. The first one is a good assumption but the second one can be 
far away from reality. Let us now present a general bottom condition formula for 
normal modes that depends on the bottom characteristics. 

The acoustic pressure due to the mode is 

p„ = ZMRnir,e)e^‘^‘ . (5.1) 

The vertical component of the momentum equation relating p„ to the vertical particle 
velocity in the mode is 



Pi 



dt 




{B.l) 



where p, is the water density and the bottom is assumed to be flat. Using (B.l) in (B.2) 
we get 



= 



-^5 

P\^ dz 



/ 



(5.3) 



At the bottom boundary, the requirement is that both p„ and are continuous 
across the water-sediment interface. In other words, the normal acoustic impedance 
across the interface is continuous. The normal acoustic impedance associated with the 
mode is 



Pn 



(5.4) 



By using (B.l) and (B.3) in (B.4) one can obtain as a general condition 



, . Pi" \ 

* ■"' .-.v, ^4-"-'’ 



(5.5) 



if at the water-sediment interface is known from the sediment properties. 
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Let us use a local plane wave approximation of the normal mode near z = II in the 
water column, i.e., 

Z„(z) = 

where A„ is a constant, R is the complex plane wave reflection coefficient and k„„ is the 
vertical wavenumber in water computed at the boundary’ (Clay and Medwin, 1977). 
With this plane wave approximatio, one can easily show that (Kinsler et al, 1982) 






Pl^lb 



( \ 2 / \2 2 

\ j \ w y 



(5.6) 



where the index 1 refers to water, 2 refers to sediment and b refers to the value at the 
boundary. R is also known as the Rayleigh Reflection Coefficient which can be equated 
as 



5 = 



Pi _ b^l2n 
P\ 

Pi ^l2n 

P\ ^2n 



2=H 



where is the vertical wavenumber in the sediment layer and is given by 



2 _ 2 

b^l2n 2 ’ 

C-) 



and K„ is the usual vertical wavenumber in the water column (Kinsler et al, 1982). 
Substituting (B.6) in (B.5) we can rewrite (B.5) as 



p< 

^2 ~Pl ^ . 



(5.7) 



A more intuitive fonn for (B.7) is 






(5.8) 



As expected, for a rigid boundaiw {R= 1), 
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dz 




= 0 



and, for a pressure release boundary' (/? = —!), 

i^r^z=H — ^ • 

B. TRAPPED MODES IN THE WATER COLUMN 

We will consider only modes that are trapped in the water layer in this WKB for- 
mulation. This means that is purely imaginary or 

f<2zn = 



with p„ real and defined by 




{B.9) 



Substituting (B.9) in (B.7) we get 

. (B.IO) 

Since (B.7) is a general boundary condition, (B.IO) is also general but is only appli- 
cable to modes trapped in the water column. We will apply (B.IO) to each class of 
normal modes in the water column. 

1. Class I 

As usual we assume that there is a turning point at z = z with 0 < z < //. In the 
exponential region (z < z < //) the solution is 

Z, = -7== cosh(0-0,) (/?. 11) 

V I f<zn I 

with 

\K^„\dz 



and 
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As (B.ll) must satisfy (B.IO), (/>(, must be 



06 = J^ \k, 



I dz— tanh 



-1 



1 d\^2n\ . Pi K 



2 1 K,^ 1 ' 



dz 



+ 



Pi 



,1 



]2=H 



(5.12) 



with /?„ given by (B.9). The only difference relative to the rigid bottom case is expression 
for All the other formulae remain unchanged. 

2. Class II 

In the oscillatory region, i < z < //, the solution is 

V '^zn '^zn 



with 



0 = 




and 



0j = I Krn • 

•'o 

(B.13) must satisfy (B.IO), implying that 



sin^ 0y/+ + tan 



-1 









~2 {D-hE)€ 



<!>, J 



= 0 



(5.14) 



where D is defined in Chapter II, Equation (22), 



Pi Pn 



l2=// • 



(B.15) 



£ = 



Pi ^:n 



Therefore, the characteristic equation following from (B.14) is 



for D <0, and 



for Z) > 0, where 





n—a 







7t— a 



a = tan 







D-\~E — (/), 
2 



All the other formulae are identical to the rigid bottom case. 

3. Class 111 

For this class the only difference relative to the rigid bottom case is that 4>t must 
now be computed with (B.12). 

4. Class IV 

The solution over the entire ocean depth is 

= — ==- sin 4> 



with 



4 > = 




To satisfy (B.IO) we require 
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for D >0, and 






dz = nn+ tan 



-I 




for D <0, with E given by (B. 15). This is now the Class IV mode characteristic equation. 
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APPENDIX C. FORTRAN PROGRAM WKBGEN 



A. PROGRAM DESCRIPTION 

The program WKBGEN is formed by a main program and several subprograms. 
The main program evaluates the minimum sound speed and then controls the subpro- 
grams. The inputs must be given in the PARAMETER statement. The inputs include 
the ocean depth (H), acoustic frequency (F), the wavenumber increment used in the k„ 
search (DK), first mode (NMI) and last one (N.MF) to be computed. The units in this 
program are those of the MKS system. The input sound speed profile is specified using 
the subprogram FUNCTION C(Z). SUBROUTINE TYPE evaluates the class of each 
mode. SUBROUTINE CHAR.AC controls the search for k„ and the turning points. The 
appropriate characteristic equation is given by FUNCTION EQCHAR and the turning 
points are evaluated by SUBROUTINE ZTURN. The auxiliary subprograms FUNC- 
TION PHASE, FUNCTION FKZ2 and FUNCTION FKZ evaluate respectively the 
phase integral, k ^„ and !k„ 1. The Airy Functions are computed by the FUNCTION'S 
A1(Z) and BI(Z) and their derivatives by FUNCTION'S DAI(Z) and DBI(Z). These four 
subprograms use the SUBROUTINE'S MMBSJR and MMBSIR in the LMSL libraries 
to evaluate the Bessel and .Modified Bessel Functions, respectively. 

As output, the mode numbers and the respective eigenvalues are printed on the 
screen. These results, as well as the depths of the turning points, are also written in a file 
whose name is specified in the CALL EXCMS statement at the begining of the main 
program. 

The numerical method used to solve for the characteristic equation and to find the 
turning depths is the simple but safe Bisection .Method (Gerald and Wheatley, 1989). 
The method to compute integrals is the Trapezoidal Rule (Gerald and Wheatley, 1989). 
Derivatives are evaluated by forw'ard or backward finite differences (Gerald and 
Wheatley, 1989). 

B. PROGRAM LISTING 

PROGRAM WKBGEN 
C 
C 

(3 ********************************** A************************************ 

c 

C MAIN PROGRAM 

C 

C 
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IMPLICIT REAL*8 (A-H.O-Z) 

PARAMETER (H= ,F= ,DK= ,NMI= ,NMF= ) 



CALL EXCMSCFILEDEF 1 DISK WKBGEN DATA A’) 

0M=8. D0*DATAN( 1. D0)*F 

ZEX=DMOD(H,10.D0) 



CS=C(0) 

CB=C(H) 



CM=99999. DO 
DO Z=0.D0,H-ZEX,10.D0 
CM=DMIN1(CM,C(Z)) 
END DO 

CM=DMIN1(CM,C(H)) 



XK0=0M/CM 

XKL=XK0 



DO N=NMI,NMF 
C 

DO XKI=XKL,0,DK 
C 

XKF=XKI+DK 

CNI=OM/XKI 

CNF=0M/XKF 

CALL TYPE(CNI,CS,CB,NTYPE) 

CALL ZTURN(CNI,H,ZEX,NTYPE,ZI1,ZI2) 

CALL ZTURN( CNF , H , ZEX , NTYPE , ZF 1 , ZF2 ) 

CALL CHARAC(OM, H, ZEX, N,XKI,XKF, NTYPE, ZI1,ZI2,ZF1,ZF2,NS0L,XKS0L, 
$ ZT1,ZT2) 

IF (NSOL. EQ. 1) GO TO 100 
C 

END DO 
C 

100 PRINT*, N,XKSOL 
C 

IF (NTYPE. EQ. 1. OR. NTYPE. EQ. 2) THEN 
WRITE (1,*) N, NTYPE, XKS0L,ZT1 
ELSE IF (NTYPE. EQ. 3) THEN 

WRITE (1,*) N, NTYPE, XKS0L,ZT1,ZT2 
ELSE 

WRITE (1,*) N, NTYPE, XKSOL 
END IF 
C 

XKL=XKSOL 
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END DO 
STOP 
END 
C 
C 
C 

C5V******5V**lV*lVlV************5V***5V5V5V*'lV*lV*****5VlV*******5V***************5V*5V* 

c 

c 

c 

SUBROUTINE TYPE(CN,CS ,CB ,NTYPE) 

REAL*8 CN,CS,CB 
C 

IF (CN. GT.CS.AND.cn. LT.CB) THEN 
NTYPE=1 

ELSE IF (CN. LT.CS.AND.cn. GT.CB) THEN 
NTYPE=2 

ELSE IF (CN.LT.CS.AND.cn. LT.CB) THEN 
NTYPE=3 
ELSE 

NTYPE=4 
END IF 
C 

RETURN 

END 

C 

C 

C 

c 

c 

c 

SUBROUTINE CHARAC ( ON , H , ZEX , NM , XI , XF , NT , Z I 1 , Z I 2 , ZF 1 , ZF2 , NSOL , XSOL , 

$ Z21 Z22) 

IMPLICIT REAL*8 (A-H,0-Z) 

XSOL=9999.DO 

NSOL=l 

FI=EQCHAR(NM,NT,H,OM,XI,ZIl,ZI2) 

FF=EQCHAR( NM , NT , H , OM , XF , ZF 1 , ZF2 ) 

IF ((FI*FF).GT. O.DO) THEN 
NS0L=0 
GO TO 100 

ELSE IF (FI. EQ. O.DO) THEN 
XSOL=XI 
GO TO 500 

ELSE IF (FF. EQ. O.DO) THEN 
XSOL=XF 
GO TO 500 
ELSE 

CONTINUE 
END IF 

DO 200 N=l,500 
XM2=(XF+XI)/2.D0 
IF(N. EQ. 1) GO TO 50 
IF ((XM2-XM1). EQ. 0. DO) THEN 
XS0L=.XM2 
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50 



GO TO 500 
END IF 
XM1=XM2 
CN2=OM/XM2 

CALL ZTURN ( CN2 , II , ZEX , NT , Z 2 1 , Z 2 2 ) 
F2=EQCHAR(NM,NT,H.OM,X2,Z21,Z22) 
IF ((F2*FF). LT. 0. DO) THEN 
XI=XM2 
ELSE 

XF=XM2 
FF=F2 
END IF 
200 CONTINUE 
100 CONTINUE 
500 CONTINUE 
RETURN 
END 






SUBROUTINE ZTURN( CN , H , ZEX , NT , ZTl , ZT2 ) 
IMPLICIT REAL*8 (A-H.O-Z) 

REAL*8 NM 

ZT1=99999. DO 

ZT2=99999. DO 

DO 100 D= 0.D0,H,10.D0 

XF=D+10. DO 

XI=D 

FI=C(XI)-CN 

FF=C(XF)-CN 

IF ((FI*FF).GT. O.DO) THEN 
GO TO 100 

ELSE IF (FI. EQ. O.DO) THEN 
XSOL=XI 
GO TO 500 

ELSE IF (FF. EQ. O.DO) THEN 
XSOL=XF 
GO TO 500 
ELSE 

CONTINUE 
END IF 

DO 200 N=l,500 
XM2=(XF+XI)/2. DO 
IF(N. EQ. 1) GO TO 50 
IF ((XM2-XM1). EQ. 0. DO) THEN 
XS0L=XM2 
GO TO 500 
END IF 
50 XM1=XM2 

F2=C(XM2)-CN 

IF ((F2*FF). LT. 0. DO) THEN 
XI=XM2 
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ELSE 

XF=XM2 
FF=F2 
EN’D IF 
200 CONTINUE 
100 CONTINUE 
500 CONTINUE 
ZT1=XS0L 

IF (NT. EQ. 3) THEN 
DO 110 D= H,0. D0,-10. DO 
XF=D-10. DO 
XI=D 

FI=C(XI)-CN 

FF=C(XF)-CN 

IF ((FI*FF).GT. O.DO) THEN 
GO TO 110 

ELSE IF (FI.EQ. O.DO) THEN 
XSOL=XI 
GO TO 510 

ELSE IF (FF.EQ. O.DO) THEN 
XSOL=XF 
GO TO 510 
ELSE 

CONTINUE 
END IF 

DO 210 N=l,510 
XM2=(XF+XI)/2.D0 
IF(N. EQ. 1) GO TO 51 
IF ((XM2-XM1). EQ. 0. DO) THEN 
XS0L=XM2 
GO TO 510 
END IF 
51 XM1=XH2 

F2=C(XM2)-CN 

IF ((F2*FF). LT. 0. DO) THEN 
XI=XM2 
ELSE 

XF=XM2 
FF=F2 
END IF 
210 CONTINUE 
110 CONTINUE 
510 CONTINUE 
ZT2=XS0L 
END IF 
RETURN 
END 



FUNCTION EQCHAR(NM,NT,H,0M,XI,ZT1,ZT2) 
IMPLICIT REAL*8 (A-H.O-Z) 
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PI=4. D0*DATAN(1. DO) 

D=(FKZ(OM,XI,H)-FKZ(OM,XI,H-5. DO) )/( 10. DO*FKZ(OM,XI ,H)**2) 
IF (DABS(D). GT. . 999D0) THEN 
D1=(D/DABS(D))*. 999D0 
ELSE 
D1=D 
END IF 



IF (NT. EQ. 1) THEN 

EQCHAR=PHASE( 0 . DO , ZTl , OM , XI ) - ( NM- . 25D0 )*PI 
IF (H-ZTl. GT. 10. DO) THEN 
PHI B=PHASE( ZTl , H , OM , XI ) -DATANH( D 1 ) 
EQCHAR=EQCHAR+DATAN(DEXP(-2. D0*PHIB)/2. DO) 

ELSE 

GAM=(FKZ2(0M,XI,H)-FKZ2(0M,XI,ZT1))/(ZT1-H) 

ZH=( GAM/ DAB S ( GAM ))*DABS( GAM )**(!. DO/3. D0)*(H-ZT1) 
A=DBI(ZH) 

B=DAI(ZH) 

EQCHAR=EQCHAR-DATAN( B/A) 

END IF 



ELSE IF (NT.EQ. 2) THEN 

EQCHAR= PHASE ( ZT 1 , H , OM , XI ) 

IFCZTl.GT. lO.DO) THEN 
PHIS=PHASE(0. DO, ZTl, OM, XI) 

ALF1=DEXP(PHIS)+D*DEXP(-1. D0*PHIS)/2. DO 
ALF2=((DEXP(-1. D0*PHIS)/2. DO) -D*DEXP(PHIS) ) 

IF (DLOG10(DABS(ALF1))-DLOG10(DABS(ALF2)). GT. 60. DO) THEN 
IF (D. GE.O) THEN 
ALF=PI/2. DO 
ELSE 

ALF=-1.D0*PI/2.D0 
END IF 
ELSE 

ALF=DATAN( ALF 1 / ALF 2 ) 

END IF 

IF (D. GE.O. DO) THEN 

EQCHAR=EQCHAR-(NM-1. 25D0)*PI+ALF 
ELSE 

EQCHAR=EQCHAR-(NM-. 25D0)*PI+ALF 
END IF 
ELSE 

GAM=(FKZ2(OM,XI,ZT1)-FKZ2(OM,XI,0. D0))/ZT1 
ZS=(GAM/DABS(GAM))*DABS(GAM)**( 1. DO/3. DO) 

A=BI(ZS) 

B=AI(ZS) 

IF (D. LT. O.DO) THEN 

EQCHAR=EQCHAR-(NM-. 25D0)*PI+DATAN( ( A+D*B)/( B-D*A) ) 

ELSE 

EQCHAR=EQCHAR-(NM-1. 25D0)*PI+DATAN((A+D*B)/(B-D*A)) 

END IF 
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END IF 



ELSE IF (NT. EQ. 3) THEN 
IF (H-ZT2.GT. 10. DO) THEN 
PHIB=PHASE( ZT2 , H , OM ,XI ) -DATANH( Dl) 
ALF2=DATAN(DEXP( -2. D0*PHIB)/2. DO) 

ELSE 

GAM=(FKZ2(0M,XI,H)-FKZ2(0M,XI,ZT2))/(ZT2-H) 

ZH=( GAM/DABSC GAM) )*DABS( GAM)**( 1. DO/3. D0)*( H-ZT2 ) 
A=DBI(ZH) 

B=-l. D0*DAI(ZH) 

ALF2=DATAN(B/A) 

END IF 

IF(ZT1. GT. 10. DO) THEN 
PHIS=PHASE(O.DO,ZT1,OM,XI) 

ALF1=DATAN(DEXP( -2. D0*PHIS)/2. DO) 

ELSE 

GAM=(FKZ2(OM,XI,ZT1)-FKZ2(OM,XI,O.DO))/ZT1 
ZS=(GAM/DABS(GAM))*DABS(GAM)**(1. do/3. D0)*ZT1 
A=BI(ZS) 

B=AI(ZS) 

ALF1=DATAN(B/A) 

END IF 

EQCHAR=PHAS E ( ZT 1 , ZT2 , OM , X I ) 

EQCHAR=EQCHAR- ( NM- . 5D0 )*PI -ALF 1+ALF2 



ELSE 

EQCHAR=PHASE( 0 . DO , H , OM , X I ) 

IF (DABS(D).GT. l.D-60) THEN 
E=1.D0/D 

IF(D. GT. 0. DO) THEN 

EQCHAR=EQCHAR-NM*PI+DATAN( E ) 
ELSE 

EQCHAR=EQCHAR-NM*P I -DATAN( E ) 
END IF 
ELSE 

EQCHAR=EQCHAR-NM*PI+PI/2. DO 
ENT) IF 



END IF 

RETURN 

END 



A*********************************************** *********************** 
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FUNCTION PHASE(A,B,OM,XKN) 

IMPLICIT REAL*8 (A-H,0-Z) 

XNL=DINT((B-A)/5. DO) 

NL=NINT(XNL) 

EX=DMOD((B-A),5.DO) 

PHASE=0. DO 
DO 1=1, NL 
XI=DBLE(I) 

PHASE=PHASE+(FKZ(OM,XKN,(XI-l. D0)*5. DO+A) 

$ +FKZ(0M,XKN,XI*5.D0+A))*2.5D0 

END DO 

PHASE=PHASE+( FKZ ( OM , XKN , XNL*5 .DO+A) 

$ +FKZ(OM,XKN,B))*EX/2.DO 

RETURN 
END 
C 
C 
C 

Q)fr*yf***)^r**Vc)^r***^V**)fr)^n^r)fr)fr)^r***)^r***)^r)^r***)fr********)fr*)^r)^r)fr*********Vf***********:fr 



c 

c 

c 



FUNCTION FKZ(OM,XKN,Z) 

REAL*8 FKZ,OM,XKN,Z 

FKZ=(OM/C(Z))**2-XKN**2 

FKZ=DSQRT( DABS( FKZ ) ) 

RETURN 

END 



C 

C 

C 

C*********************************************************************** 

c 

c 

c 

FUNCTION FKZ2(OM,XKN,Z) 

REAL*8 FKZ2,OM,XKN,Z 
FKZ2=(OM/C(Z))**2-XKN**2 
RETURN 
END 
C 
C 
C 

C*yf*yr*yf**yryfyfyc**yr****yr*************yc*******yr****yr***yr******yryp*yrycyr*****yr*yr 

c 

c 

c 



FUNCTION C(Z) 
REAL*8 C,Z 
C= 

RETURN 

END 



C 

C 

C 

Qyfyfyc*yryr***yryryryr***yf*vc****vcvc*******vcyr**vf*yryr*******vc****vcyryfyf*yryfyfyrvcyf***yc**** 
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FUNCTION AI(Z) 

IMPLICIT REAL*8 (A-H.O-Z) 
REAL*8 J13,JM13,I13,IM13 
DIMENSION RJ(2),WK(4),B(2) 



IF (Z.EQ. O.DO) THEN 
AI=. 35502805D0 
RETURN 

ELSE IF (Z.LT. O.DO) THEN 



ARG=2. D0*(DSQRT( -1. D0*Z)**3)/3. DO 
N=2 



0RDER=2. DO/ 3. DO 

CALL MMBS JR( ARC , ORDER , N , RJ , WK , lER) 
JM13=4. D0*RJ( 1 ) / ( 3 . D0*ARG) -RJ( 2 ) 



0RDER=1. DO/ 3. DO 

CALL MMBS JR( ARG , ORDER , N , R J , WK , lER ) 
J13=RJ(1) 



AI=DSQRT(-1. D0*Z)*( JM13+J13)/3. DO 



RETURN 

ELSE 



ARG=2. D0*(DSQRT(Z)**3)/3. DO 

0RDER=2. DO/ 3. DO 

NB=2 

I0PT=1 



CALL MMBS IR( ARG , ORDER , NB , lOPT , B , lER) 
IM13=4. D0*B( 1) /( 3. D0*ARG)+B( 2) 



0RDER=1. DO/3. DO 
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CALL MHBS IR( ARG , ORDER , NB , lOPT , B , I ER ) 
I13=B(1) 



AI=DSQRT(Z)*(IM13-I13)/3. DO 
RETURN 



END IF 



END 



FUNCTION BI(Z) 

IMPLICIT REAL*8 (A-H.O-Z) 
REAL*8 J13,JM13,I13,IM13 
DIMENSION RJ(2),WK(4),B(2) 



IF (Z.EQ. O.DO) THEN 
BI=. 61492663D0 
RETURN 

ELSE IF (Z. LT. O.DO) THEN 



ARG=2, D0*(DSQRT(-1. D0*Z)**3)/3. DO 
N=2 



0RDER=2.D0/3.D0 

CALL MMBSJRC ARG, ORDER, N,RJ,WK,IER) 
JM13=4. D0*RJ( l)/(3. D0*ARG)-RJ(2) 



0RDER=1. DO/3. DO 

CALL MMBSJRC ARG , ORDER , N , R J , WK , I ER) 
J13=RJ(1) 



BI=DSQRT( -1. D0*Z/3. D0)*( JM13-J13) 
C 
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RETURN 

ELSE 



ARG=2. D0*(DSQRT(Z)**3)/3. DO 

0RDER=2. DO/3. DO 

NB=2 

I0PT=1 



CALL MMBS IR( ARC , ORDER , NB , lOPT , B , lER) 
IM13=4. D0*B( l)/(3. D0*ARG)+B(2) 



0RDER=1. DO/3. DO 

CALL MMBS IR( ARG , ORDER , NB , lOPT , B , lER) 
I13=B(1) 



BI=DSQRT(Z/3.D0)*(IM13+I13) 

RETURN 



END IF 



END 



FUNCTION DAI(Z) 

IMPLICIT REAL*8 (A-H,0-Z) 
REAL*8 J23,JM23,I23,IM23 
DIMENSION RJ(2),WK(4),B(2) 



IF (Z.EQ. O.DO) THEN 
DAI=-. 25881940D0 
RETURN 

ELSE IF (Z.LT. O.DO) THEN 
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ARG=2. D0*(DSQRT(-1. D0*Z)**3)/3. DO 
N=2 



0RDER=1. DO/3. DO 

CALL MMBS JR( ARC , ORDER , N , R J , WK , lER ) 
JM23=2. D0*RJ( l)/(3. D0*ARG)-RJ(2) 



0RDER=2. DO/3. DO 

CALL MMBS JR( ARG , ORDER , N , RJ , WK , lER) 
J23=RJ(1) 



DAI=Z*(JM23-J23)/3. DO 



RETURN 

ELSE 



ARG=2. D0*(DSQRT(Z)**3)/3. DO 

0RDER=1. DO/3. DO 

NB=2 

I0PT=1 



CALL MMBS IR( ARG , ORDER ,NB , lOPT, B , lER) 
IM23=2. D0*B( 1) /( 3. D0*ARG)+B( 2) 



0RDER=2. DO/3. DO 

CALL MMBS IR( ARG , ORDER , NB , lOPT , B , lER ) 
I23=B(1) 



DAI=-1. D0*Z*( IM23-I23)/3. DO 
RETURN 



END IF 



END 
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c 



FUNCTION DBI(Z) 

IMPLICIT REAL*8 (A-H.O-Z) 
REAL*8 J23,JM23,I23,IM23 
DIMENSION RJ(2) ,WK(4),B(2) 



IF (Z.EQ. O.DO) THEN 
DBI=. 44828836D0 
RETURN 

ELSE IF (Z.LT. 0. DO) THEN 



ARG=2. D0*(DSQRT( -1. D0*Z)**3)/3. DO 
N=2 



ORDER=1.DO/3.DO 

CALL MMBS JR( ARC , ORDER , N , R J , WK , lER) 
JM23=2. D0*RJ( 1 ) /( 3. DO* ARC) -RJ( 2 ) 



ORDER=2.DO/3.DO 

CALL MMBS JR( ARC , ORDER , N , RJ , WK , lER) 
J23=RJ(1) 



DBI=-1. D0*Z*(JM23+J23)/DSQRT(3. DO) 



RETURN 

ELSE 



ARG=2. D0*(DSQRT(Z)**3)/3.D0 

0RDER=1. DO/3. DO 

NB=2 

I0PT=1 



CALL MMBSIRC ARG , ORDER , NB , lOPT, B , lER) 
IM23=2. D0*B( l)/(3. D0*ARG)+B(2) 
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c 

c 



0RDER=2. DO/3. DO 

CALL M.MBS IR( ARC , ORDER , NB , lOPT , B , lER) 
I23=B(1) 



DBI=Z*( IM23+I23)/DSQRT( 3. DO) 
RETURN 

END IF 

END 
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